% GAUTIER LE BIHAN - 2020
% Replication files for "Shocks vs Menu Costs: Patterns of Price Rigidity in an Estimated Multi-Sector Menu-Cost Model?" Review of Economics and Statistics
%
% This code produces  Figure K + Table B

clear all;
tic

cd ..\..\Simulations_VC\MS_produits_CPlus

load ..\..\Simulations_VC\MS_produits_CPlus\stat_d
param=stat_d;

load ..\..\Estimation_param\MS_produits_CPlus\actual_moments_k
weight=actual_moments_k(:,12);
sum(weight);

secteur=actual_moments_k(:,1);

f=param(:,5);
%muc=param(:,2);
f_p=param(:,6);
med=param(:,7);
interq=param(:,8);
kur=param(:,9);
q1=param(:,10);
q2=param(:,11);
q3=param(:,12);

moments_sim=[f f_p med interq kur];
moments_actual=[actual_moments_k(:,2:6)];
size_param=size(moments_sim,2);

fprintf("Mean Simulated Moments across Products")
mom_mean=(weight'*moments_sim)/sum(weight)
kf_sim=(moments_sim(:,5)./moments_sim(:,1));
kf_act=(moments_actual(:,5)./moments_actual(:,1));
kf_sims=(weight'*kf_sim)/sum(weight)
kf_acts=(weight'*kf_act)/sum(weight)

mom_med=ones(1,size_param+2);
stat_sect=moments_sim;
%figure;hist(param_sect(:,3));
v=size(stat_sect,1);
weight_sect=weight;
 if v>0
    for i=1:size(stat_sect,2)
        [sortx,order] = sort(stat_sect(:,i));
        sortw = weight_sect(order);

         midpoint = sum(sortw)/2;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              mom_med(i) = mean(sortx([j j+1]));
            else
                
                mom_med(i) = sortx(j+1);
            end
     mom_med(end)= csumw(end);
     mom_med(end-1)= v;
    end  
 end
fprintf("Median Simulated Moments across Products")
mom_med


figure(1);
%P0
subplot(2,3,1);
x=moments_sim(:,1);
y=moments_actual(:,1);


scatter(y, x);

x1=[0  0.2];
y1=x1;
line(x1,y1);

axis([x1 x1]);
ylabel('Freq. sim.') % x-axis label
xlabel('Freq. data') % y-axis label


%MU_c
subplot(2,3,2);
x=moments_sim(:,2);

%x=x(x>0.4);
y=moments_actual(:,2);

%y=y(y>0.2);

scatter(y, x);
x1=[0.4  1];
y1=x1;
line(x1,y1);

axis([x1 x1]);
ylabel('Frac inc. sim.') % x-axis label
xlabel('Frac inc. data') % y-axis label

%MU_c
subplot(2,3,4);
x=moments_sim(:,3);

%x=x(x>0.4);
y=moments_actual(:,3);

%y=y(y>0.2);

scatter(y, x);
x1=[0 0.15];
y1=x1;
line(x1,y1)

axis([x1 x1]);
ylabel('Median sim.') % x-axis label
xlabel('Median data') % y-axis label

%MU_c
subplot(2,3,5);
x=moments_sim(:,4);

%x=x(x>0.4);
y=moments_actual(:,4);

%y=y(y>0.2);

scatter(y, x);
x1=[0 0.25];
y1=x1;
line(x1,y1);
axis([x1 x1])
ylabel('Interq. sim.') % x-axis label
xlabel('Interq. data') % y-axis label

%MU_c
subplot(2,3,6);
x=moments_sim(moments_sim(:,5)>0 & moments_sim(:,5)<10 &  moments_actual(:,5)<10,5);

%x=x(x>0.4);
y=moments_actual(moments_sim(:,5)>0 & moments_sim(:,5)<10 &  moments_actual(:,5)<10,5);

%y=y(y>0.2);
axis([0 10 0 8]);
scatter(y, x);
x1=[0 6];

y1=x1;
line(x1,y1);

lm = fitlm(y,x)
%p(i,:)=lm.Coefficients(2,1)
%v(i,:)=lm.Coefficients(2,4)
x2=[0 8];

ylabel('Kurtosis sim.') % x-axis label
xlabel('Kurtosis data') % y-axis label

 print('..\..\..\figures\scatter_act_sim.pdf','-dpdf', '-fillpage')
 print('..\..\..\figures\scatter_act_sim','-depsc')






stat_med_sect=ones(6,size_param+2);
for k=1:6
stat_sect=moments_sim((secteur==k), :);
%figure;hist(param_sect(:,3));
v=size(stat_sect,1);
weight_sect=weight((secteur==k),:);
 if v>0
    for i=1:size(stat_sect,2)
        [sortx,order] = sort(stat_sect(:,i));
        sortw = weight_sect(order);

         midpoint = sum(sortw)/2;
         csumw = cumsum(sortw);
         j = find(csumw<=midpoint,1,'last');
         
        if isempty(j)==1
             j=1;
         end
            if csumw(j) == midpoint
              stat_med_sect(k,i) = mean(sortx([j j+1]));
            else
                
                stat_med_sect(k,i) = sortx(j+1);
            end
     stat_med_sect(k,end)= csumw(end);
     stat_med_sect(k,end-1)= v;
    end  
 end
end
fprintf("By sector");
stat_med_sect
